home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ASTRNOMY / AA_51.ZIP / OJUPITER.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  5KB  |  240 lines
  1. /* Orbital elements and perturbations for the planet Jupiter
  2.  * using formulas given by Meeus
  3.  */
  4. #include "planet.h"
  5. #include "kep.h"
  6.  
  7.  
  8. /* Calculations needed for Jupiter through Neptune:
  9.  */
  10. int pjup()
  11. {
  12.  
  13. nu = T/5.0 + 0.1;
  14. P = 237.47555 + 3034.9061*T;
  15. Q = 265.91650 + 1222.1139*T;
  16. S = 243.51721 + 428.4677*T;
  17. f = 5.0*Q - 2.0*P;
  18. V = mod360(f)*DTR;
  19. f = 2.0*P - 6.0*Q + 3.0*S;
  20. W = mod360(f)*DTR;
  21. f = Q - P;
  22. ze = mod360(f)*DTR;
  23. f = S - Q;
  24. psi = mod360(f)*DTR;
  25. P = mod360(P)*DTR;
  26. Q = mod360(Q)*DTR;
  27. S = mod360(S)*DTR;
  28.  
  29. sinV = sin(V);        
  30. cosV = cos(V);
  31. /* Use trigonometry to speed up sin, cos of multiple angles.
  32.  * This reduces the relative accuracy, but we only need
  33.  * absolute accuracy to about 8 places.
  34.  */
  35. sin2V = 2.0*sinV*cosV;
  36. cos2V = cosV*cosV - sinV*sinV;
  37. sinW = sin(W);
  38. cosze = cos(ze);
  39. sinze = sin(ze);
  40.  
  41. sin2ze = 2.0*sinze*cosze;
  42. cos2ze = cosze*cosze - sinze*sinze;
  43.  
  44. sin3ze = sinze*cos2ze + cosze*sin2ze;
  45. cos3ze = cosze*cos2ze - sinze*sin2ze;
  46.  
  47. sin4ze = sinze*cos3ze + cosze*sin3ze;
  48. cos4ze = cosze*cos3ze - sinze*sin3ze;
  49.  
  50. sin5ze = sinze*cos4ze + cosze*sin4ze;
  51. cos5ze = cosze*cos4ze - sinze*sin4ze;
  52.  
  53. sinQ = sin(Q);
  54. cosQ = cos(Q);
  55.  
  56. sin2Q = 2.0*sinQ*cosQ;
  57. cos2Q = cosQ*cosQ - sinQ*sinQ;
  58.  
  59. sin3Q = sinQ*cos2Q + cosQ*sin2Q;
  60. cos3Q = cosQ*cos2Q - sinQ*sin2Q;
  61. return(0);
  62. }
  63.  
  64. int djup();
  65.  
  66. /* Elements and perturbations for Jupiter
  67.  */
  68. int ojupiter(e,J)
  69. struct orbit *e;
  70. double J;
  71. {
  72.  
  73. e->epoch = J;
  74. manoms(J);
  75. pjup();
  76. e->M = M5;
  77. e->a = 5.202561;
  78. e->ecc = (( -0.0000000017*T - 0.0000004676)*T + 0.000164180)*T + 0.04833475;
  79. #if OFDATE
  80. e->equinox = J;
  81. e->i = ( 0.0000039*T - 0.0056961)*T + 1.308736;
  82. f = (( 0.00000508*T + 0.00070405)*T + 0.5994317)*T + 273.277558;
  83. e->w = mod360(f);
  84. f = (( -0.00000851*T + 0.00035222)*T + 1.0105300)*T + 99.443414;
  85. e->W = mod360(f);
  86. f = (( -0.00000165*T + 0.0003347)*T + 3036.301986)*T + 238.049257;
  87. e->L = mod360(f);
  88. #else
  89. e->equinox = J2000;
  90. e->i = ((1.27e-7*T + 2.942e-5)*T - 0.0022374)*T + 1.305288;
  91. f = ((8.999e-6*T - 0.00021857)*T + 0.0478404)*T + 273.829584;
  92. e->w = mod360(f);
  93. f = (( -1.2460e-5*T + 0.00096672)*T + 0.1659357)*T + 100.287838;
  94. e->W = mod360(f);
  95. f = e->M + e->w + e->W;
  96. e->L = mod360(f);
  97. #endif
  98. djup(e);
  99. return(0);
  100. }
  101.  
  102.  
  103.  
  104. /* Program split arbitrarily to get it through
  105.  * Microsoft C compiler:
  106.  */
  107. int djup(e)
  108. struct orbit *e;
  109. {
  110. double p, q;
  111.  
  112. /* perturbations in mean longitude */
  113. p = ((-0.004692*nu - 0.010281)*nu + 0.331364)*sinV
  114.    +((0.002075*nu - 0.064436)*nu + 0.003228)*cosV
  115.    -(( -0.000489*nu + 0.000275)*nu + 0.003083)*sin2V
  116.    +0.002472*sinW
  117.    +0.013619*sinze
  118.    +0.018472*sin2ze
  119.    +0.006717*sin3ze
  120.    +0.002775*sin4ze;
  121. f = (0.007275 - 0.001253*nu)*sinze
  122.     +0.006417*sin2ze
  123.     +0.002439*sin3ze
  124.     -(0.033839 + 0.001125*nu)*cosze
  125.     -0.003767*cos2ze;
  126. p += f*sinQ;
  127. f = -(0.035681 + 0.001208*nu)*sinze
  128.      -0.004261*sin2ze
  129.      +0.002178
  130.      +(-0.006333 + 0.001161*nu)*cosze
  131.      -0.006675*cos2ze
  132.      -0.002664*cos3ze;
  133. p += f*cosQ;   
  134. f =  -0.002572*sinze
  135.      -0.003567*sin2ze;
  136. p += f*sin2Q;
  137. f =  0.002094*cosze
  138.     +0.003342*cos2ze;
  139. p += f*cos2Q;
  140.  
  141. /* perturbations in the perihelion */
  142. q = (0.007192 - 0.003147*nu)*sinV
  143.    +((0.000197*nu - 0.000675)*nu - 0.020428)*cosV;
  144. f = (0.007269 + 0.000672*nu)*sinze
  145.     -0.004344
  146.     +0.034036*cosze
  147.     +0.005614*cos2ze
  148.     +0.002964*cos3ze;
  149. q += f*sinQ;
  150. f =  0.037761*sinze
  151.     +0.006158*sin2ze
  152.     -0.006603*cosze;
  153. q += f*cosQ;
  154. f = -0.005356*sinze
  155.     +0.002722*sin2ze
  156.     +0.004483*cosze
  157.     -0.002642*cos2ze;
  158. q += f*sin2Q;
  159. f =  0.004403*sinze
  160.     -0.002536*sin2ze
  161.     +0.005547*cosze
  162.     -0.002689*cos2ze;
  163. q += f*cos2Q;        
  164. e->w += q;
  165. #if DEBUG
  166. printf( "pperih %.5e\n", q );
  167. #endif
  168. f = p - (q/(e->ecc));
  169. e->M += f;
  170. #if DEBUG
  171. printf( "pmanom %.5e\n", f );
  172. #endif
  173. e->L += p;
  174. #if DEBUG
  175. printf( "plong %.5e\n", p );
  176. #endif
  177.  
  178. /* perturbations in the eccentricity */
  179. q = ((-.0000043*nu + .0000130)*nu + .0003606)*sinV;
  180. q += (.0001289 - .0000580*nu)*cosV;
  181. f = -0.0006764*sinze
  182.     -0.0001110*sin2ze
  183.     -0.0000224*sin3ze
  184.     -0.0000204
  185.     +( 0.0001284 + 0.0000116*nu)*cosze
  186.     +0.0000188*cos2ze;
  187. q += f*sinQ;
  188. f = (0.0001460 + 0.0000130*nu)*sinze
  189.     +0.0000224*sin2ze
  190.     -0.0000817
  191.     +0.0006074*cosze
  192.     +0.0000992*cos2ze
  193.     +0.0000508*cos3ze
  194.     +0.0000230*cos4ze
  195.     +0.0000108*cos5ze;
  196. q += f*cosQ;
  197. f = -(0.0000956 + 0.0000073*nu)*sinze
  198.     +0.0000448*sin2ze
  199.     +0.0000137*sin3ze
  200.     +( -0.0000997 + 0.0000108*nu)*cosze
  201.     +0.0000480*cos2ze
  202.     +0.0000148*cos3ze;
  203. q += f*sin2Q;
  204. f = ( -0.0000956 +0.0000099*nu)*sinze
  205.     +0.0000490*sin2ze
  206.     +0.0000158*sin3ze
  207.     +0.0000179
  208.     +( 0.0001024 + 0.0000075*nu)*cosze
  209.     -0.0000437*cos2ze
  210.     -0.0000132*cos3ze;
  211. q += f*cos2Q;     
  212. e->ecc += q;
  213.  
  214. /* perturbations in the semimajor axis */
  215. q = -0.000263*cosV
  216.     +0.000205*cosze
  217.     +0.000693*cos2ze
  218.     +0.000312*cos3ze
  219.     +0.000147*cos4ze;
  220. f =  0.000299*sinze
  221.     +0.000181*cos2ze;
  222. q += f*sinQ;
  223. f =  0.000204*sin2ze
  224.     +0.000111*sin3ze
  225.     -0.000337*cosze
  226.     -0.000111*cos2ze;
  227. q += f*cosQ;        
  228. e->a += q;
  229. return(0);
  230. }
  231.  
  232. int cjupiter(e)
  233. struct orbit *e;
  234. {
  235.  
  236. e->plat = 0.0;
  237. return(0);
  238. }
  239.  
  240.